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Abstract 

Background: The computational methods provide condition for investigation related to the process of drug 
delivery, such as convection and diffusion of drug in extracellular matrices, drug extravasation from microvessels or 
to lymphatic vessels. The information of this process clarifies the mechanisms of drug delivery from the injection 
site to absorption by a solid tumor. In this study, an advanced numerical method is used to solve fluid flow and 
solute transport equations simultaneously to investigate the effect of tumor shape and size on drug delivery to 
solid tumor. 

Methods: The advanced mathematical model used in our previous work is further developed by adding solute 
transport equation to the governing equations. After applying appropriate boundary and initial conditions on 
tumor and surrounding tissue geometry, the element-based finite volume method is used for solving governing 
equations of drug delivery in solid tumor. Also, the effects of size and shape of tumor and some of tissue transport 
parameters such as effective pressure and hydraulic conductivity on interstitial fluid flow and drug delivery are 
investigated. 

Results: Sensitivity analysis shows that drug delivery in prolate shape is significantly better than other tumor 
shapes. Considering size effect, increasing tumor size decreases drug concentration in interstitial fluid. This study 
shows that dependency of drug concentration in interstitial fluid to osmotic and intravascular pressure is negligible. 

Conclusions: This study shows that among diffusion and convection mechanisms of drug transport, diffusion is 
dominant in most different tumor shapes and sizes. In tumors in which the convection has considerable effect, the 
drug concentration is larger than that of other tumors at the same time post injection. 
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Background 

Cancer, the main cause of death in developed countries, 
is the second leading cause of death in developing coun- 
tries [1]. Solid tumors account for 85% of human can- 
cers [2]. Chemotherapy is one of the ways widely used 
for cancer treatment. Based on the findings from clinical 
applications, most cancer treatments with drugs fail to 
eliminate solid tumors completely [3]. The computational 
method can investigate why systemic administration 
cannot distribute drug uniformly in tumors. The drug 
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exchange between microvessels and extracellular matrices, 
drug removal by lymphatic system, drug diffusion and 
convective transport in extracellular matrices should be 
included by mathematical simulation. Computational fluid 
dynamics (CFD) can model the whole drug delivery 
process and clarify the mechanisms of drug delivery from 
the injection site to absorption by a solid tumor. 

Baxter and Jain, based on the theoretical framework 
in their ID mathematical method, found the effective 
factors on drug delivery such as microvessel perme- 
ability, interstitial fluid pressure (IFP), and interstitial 
fluid velocity (IFV) [4-7]. 

Improving the ID model of Baxter et al. [5,6,8] and 
Saltzman et al. [9], Wang et al. [10-12] developed a 
simulation framework of drug delivery to tumors by 



o 



BioMed Central 



© 2014 Sefidgar et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
Commons Attribution License (http://creativecommons.Org/licenses/by/4.0), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain 
Dedication waiver (http://creativecommons.0rg/publicdomain/zero/l.O/) applies to the data made available in this article, 
unless otherwise stated. 



Sefidgar et al. Journal of Biological Engineering 2014, 8:12 
http://www.jbioleng.org/content/871 /1 2 



Page 2 of 13 



considering the complex 3D geometry. Wang and Li 
[10] used modified MRI images for tumor geometry. 
They considered interstitial fluid flow with blood and 
lymphatic drainage in their model. Wang et al. [11] 
studied the effect of elevated interstitial pressure, con- 
vective flux, and blood drainage on the delivery of speci- 
fied solute to brain tumors. 

The study of tissue transport property effect on drug 
delivery is considered in recent studies. Zhao et al. [13] 
used a 3D computational model to predict the distribu- 
tion of IFV, IFP, and solute transport through a tumor. 
Arifin et al. [14] studied the sensitivity of drug distribu- 
tion to physiochemical properties in realistic models of 
brain tumors. A specific tumor captured by MRI is used 
by Pishko et al. [15] for modeling drug distribution in tis- 
sue with spatially-varying porosity and vascular permeabil- 
ity. The sensitivity of solute distribution to tumor shape 
and size is not considered in above mentioned works. 

In our previous work [2], tumor shape and size effect 
on drug delivery is investigated by modeling interstitial 
fluid flow and assuming that drug particles flow with the 
interstitial fluid. In the present work, by adding the sol- 
ute transport equation to the previous developed model 
in our group [16-20], new governing equations are in- 
vestigated to find drug concentration in interstitial fluid 
(DCIF). Solving fluid flow and solute transport equa- 
tions simultaneously, the effects of tumor shape, size, 
and tissue transport properties on drug delivery to solid 
tumor are also investigated. 

Spherical and non-spherical tumors and their surroun- 
ding normal tissue are modeled with assumption of rigid 
porous media. The vasculature as a source term and 
lymphatic vessel as a sink term vary spatially. In the fol- 
lowing parts of this paper, the sensitivity analysis pro- 
vides a better understanding of the effects of tissue 
transport parameters on drug delivery. 

Results 

Simulation of interstitial fluid flow for baseline values 
(Table 1) predicts that IFP has its greatest value in the 
tumor center. IFP is non-dimensionalized by effective 
pressure. The effective pressure, P e p is a parameter de- 
fined by intravascular pressure, plasma osmotic pressure, 
and interstitial osmotic pressure by Equation (1). The 
non-dimensionalized pressure is defined by Equation (2) 



Table 1 Interstitial transport properties used in numerical 
simulations 



Pi 



P L = 

Peff Pi-a^nl-nf) 



(1) 



(2) 



Parameter 

Lp[cm/mmHg s] 

Normal 

Tumor 
K[cm 2 /mmHg s] 

Normal 

Tumor 
S/V[cm" 1 ] 

Normal 

Tumor 
P B [mmHg] 

Normal 

Tumor 
TT B [mmHg] 

Normal 

Tumor 
7Tj[mmHg] 

Normal 

Tumor 

o 

Normal 
Tumor 

Pl 

Normal 
L PL S L /V [1/mmHg s] 
Normal 



Baseline value 



0.36X10" 7 
2.80 x 10" 7 

8.53 X10" 9 
4.13 x 10" 8 



70 
200 

15.6 
15.6 

20 
20 

10 
15 

0.91 
0.82 



1.33x10" 



Reference 

[16] 

[16] 

[16] 

[16] 

[16] 

[16] 

[16] 

[15] 
[15] 



where t stands for tumor tissue. Parameters used in 
Equation (1) and (2) are introduced in the method 
section. 



Non-dimensionalized IFP along transverse and vertical 
lines (shown in Figure 1), are illustrated in Figures 2 and 
3. The maximum value of IFP occurs in the tumor cen- 
ter (Figure 2). This value is equal to P e ff, 1.53 kPa, except 
for the case of A = 0.1 (Prolate shape). IFP has uniform 
distribution at tumor region and in the inner boundary 
(for more detail see Figure 1 and boundary condition 
section) IFP falls down sharply to around 150 Pa as 
shown in Figure 2 for all shapes except for prolate one 
(A = 0.1). In normal tissue, pressure has uniform distri- 
bution close to outer boundary (for more detail, see 
Figure 1 and boundary condition section) then it decreases 
smoothly to peripheral pressure. For prolate shape as 
shown in Figure 2, IFP reduces smoothly from tumor 
center to the outer boundary, and IFP has a non-zero 
gradient in whole domain. Figure 2 shows the IFP along 
vertical line. Results are very similar to IFP results in 
Figure 2. Only for a prolate shape IFP has different pat- 
tern from what is observed in Figure 2. Maximum value 
of IFP in the prolate shape does not occur in the tumor 
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Dirichlet boundary- 
condition for IFP. 



Inner boundary 



Continuity 
boundary condition 




Axisymmetric 
boundary condition 



Transverse line 



Outer boundary 



Figure 1 Schematic of considered geometry and boundary 
conditions. The transverse and vertical lines are used to show 
results along them. 



center and takes place somewhere between the tumor 
center and tumor periphery. 

IFV distribution along the vertical and transverse lines 
is presented in Figures 4 and 5. Maximum value of IFV 
occurs close to the inner boundary. Also, in normal tis- 
sue, IFV reaches zero far from the inner boundary. How- 
ever, for the prolate shape, velocity is not zero, especially 
along transvers line (Figure 4). 

DCIF is simulated in two cases of injection. In the first 
case, the continuous injection which leads to constant 
plasma concentration (C p = constant) is considered and 
in the second case, the bolus injection in which the 
plasma concentration decreases with time exponentially 
(C p = c£ e - ln (o.5)*/r) is stu died, in which r is the drug 
half-life in plasma (Table 2). DCIF are non-dimensionalized 
by C P for continuous injection and C p for bolus injection, 
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Figure 3 Interstitial fluid pressure profile along vertical line. 

The 1 in "x" axis shows the boundary of tumor (inner boundary). 



respectively. Average of non-dimensionalized DCIF for 
two injection cases in different times is shown in Figures 6 
and 7. DCIF of prolate shape (A = 0.1) has the maximum 
value and DCIF of oblate shape (A = 10) has the minimum 
value. Other considered shapes show the similar transient 
behavior. 

Non-dimensionalized DCIF along two lines (transverse 
and vertical) is illustrated in Figures 8 and 9. The bolus 
injection results are presented in 8 hr post injection in 
which the concentration is maximum based on Figure 7 
and for continuous injection is presented at final time of 
simulation, in 72 hr post injection. The profiles of two 
types of injections are similar in spite of different values of 
DCIF. As observed in IFV and IFP profiles, the prolate 
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Figure 2 Interstitial fluid pressure profile along transverse line. 

The 1 in "x" axis shows the boundary of tumor (inner boundary). 
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Figure 4 Interstitial fluid velocity profile along transverse line. 

The 1 in "x" axis shows the boundary of tumor (inner boundary). 
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Non-dimensionalized radius 

Figure 5 Interstitial fluid velocity profile along vertical line. The 

1 in "x" axis shows the boundary of tumor (inner boundary). 



V J 



shape is different from the other shapes in DCIF distribu- 
tion. A bump is observed in DCIF curves at the inner 
boundary of all tumor shapes. The normal tissue has uni- 
form distribution of DCIF except near the boundaries 
(inner and outer). DCIF distribution in normal tissue is 
the same for all tumor shapes. 

2D contours of DCIF in tumor region for bolus injection 
in 8 hr post injection are shown in Figure 10 for all tumor 
shapes. Results show that DCIF for tumors with A = 0.1 and 
A = 10 has different distributions. Generally, in the inner 
boundary, DCIF has its maximum value. In Figures 11 and 
12, Peclet number distribution along two lines for con- 
tinuous injection is shown. Peclet number demonstrates 
the ratio of convection to diffusion across the microvessel 



Table 2 Parameters of solute transport model used in 
numerical simulations 



Parameter 


Baseline value 


Reference 


Of 






Normal 


0.9 


[8] 


Tumor 


0.9 




Deff [m 2 /s] 






Normal 


0.1 6x 10" 12 


[8] 


Tumor 


2.0x 10" 12 




P[m/s] 






Normal 


2.2x 10~ 10 * 


[8] 


Tumor 


17.3X 10~ 10 * 




t [hrf 






Plasma 


6.1 


[8] 



*- the microvessel permeability coefficient is 10% of effective permeability 
introduced in [8]. 

**- the drug half-life in plasma introduced in results section. 




Time (hr) 

Figure 6 Average non-dimensionalized DCIF in tumor region in 

during time for continuous injection. 

I J 

wall. Results show that in the tumor region Peclet 
number is zero except for prolate shape. Peclet number 
for the prolate shape especially along short radius 
(transverse line) is greater than zero. 

Some of tissue transport parameters mentioned (ef- 
fective pressure, hydraulic conductivity, and tumor size) 
in Table 1 are investigated for sensitivity analysis. The 
values of these parameters are selected near the ranges 
reported in the literature [2,15,21]. 

Figure 13 shows the influence of changing effective 
pressure, P e p on numerical results of tumor shapes with 
the same volume. Maximum value of IFP in tumor region 
(Figure 13a) linearly increases when P e ff increases. Only in 




Time (hr) 

Figure 7 Average non-dimensionalized DCIF in tumor region in 
during time for bolus injection. 
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Figure 8 Non-dimensionalized DCIF profile along transverse line for normal. Left for bolus injection and right for continuous injection. The 
1 in "x" axis shows the boundary of tumor (inner boundary). 



the prolate shape the pressure increases more than other 
tumor shapes. Average of IFV on inner boundary also 
has the same pattern and linearly changes with P e ff 
(Figure 13b). The changes in average of DCIF in tumor 
region for two cases of injection are shown in Figure 13c 
and d. For these two cases of injection, the prolate shape 
is more sensitive to P e ff changes than other tumor shapes 
and DCIF changes around 20% when the P e ff changes from 
500 Pa to 2500 Pa. The other shapes do not show signifi- 
cant variations in DCIF for these ranges of P e ff. 

Different tumor sizes are studied with changes in their 
volume. One of important metric of disease development 
and response to tumor therapy with drug is volume of 
tumors [22-25]. To quantify response to therapeutic reg- 
imens and also assess disease progression, tumor volume is 
used as a metric in many studies, such as Char et al. [26], 
Jensen et al. [27] and Gass et al. [28]. 

As shown in Figure 14a, IFP has less value than P e ff 
when the tumor volume is smaller than 1 cm 3 . When the 
tumor volume is in the order of 1 cm 3 , IFP reaches P e ff 



and by increasing the tumor radius, IFP remains constant. 
Average of IFV on the inner boundary generally decreases 
with increasing tumor size, Figure 14b. Only in prolate 
shape in small radii, IFV increases with the tumor size. As 
shown in Figure 14c and d, the mean values of DCIF have 
the greatest value in the smallest tumor. Also, the prolate 
shape has the highest value of DCIF in all studied tumors. 

Figure 15 shows the sensitivity of IFF parameters and 
DCIF to hydraulic conductivity changes. Results show 
that in all tumor shapes, if hydraulic conductivity in- 
creases, the maximum value of IFP remains constant 
and then decreases sharply. Average of IFV increases by 
increasing hydraulic conductivity. DCIF increases by in- 
creasing hydraulic conductivity and then reaches a con- 
stant value in spite of increasing hydraulic conductivity. 

Discussion 

This study presents DCIF, IFP, and IFV in solid tumors 
surrounded by normal tissue in two types of injection; 
bolus and continuous one. The model used in this study 




Non-dimensionalized radius Non-dimensionalized radius 

Figure 9 Non-dimensionalized DCIF profile along vertical line. Left for bolus injection and right for continuous injection. The 1 in "x" axis 

shows the boundary of tumor (inner boundary). 
\ J 
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investigates the effect of two characteristics of tumors 
on concentration distribution; tumor shape and size. 

Results of high IFP in tumors are discussed in our previ- 
ous studies [2,16,17] and in the experimental results of Ari- 
fin et al [29] and Huber et al [30]. Maximum value of IFP 
in spherical tumors (1 cm radius) for baseline values in 
Table 1 is 1529.5 Pa which is close to the studies of Jain et al. 
[31], Chauhan et al. [32], and Arfin et al. [33]. The current 
results are verified by experimental data of IFP measured 
by Nielsen et al. [34]. In their work, the wick-in-needle 
technique is used to measure IFP in two types of tumors 
with the same volume as tumors in the current study. They 
found IFP in the range of 1400 Pa to 1600 Pa. 

IFV on tumor boundary in spherical tumors with base- 
line values is around 0.05 (Amis which is at the same order 



of the prediction of Jain et al. [31] and experimental obser- 
vation of Hompland et al. [35]. Also, the profile of drug 
concentration for simulation with baseline values for 
spherical tumor in different times (Figures 6 and 7) agrees 
well with Baxter and Jains predictions [8]. Results show 
that DCIF for prolate shape (A = 0.1) always has the great- 
est value. Results of transient DCIF (Figures 6 and 7) show 
that drug delivery is much easier in the prolate tumors. 
Also, the oblate shape (A = 10) has the most resistance to 
drug delivery. Following of this section, the reason of these 
phenomena is investigated. 

The uniform values of IFP for all tumor shapes except 
prolate one (Figures 2 and 3) is equal to P e ff. Non-uniform 
IFP in prolate shape results in the maximum value of 
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DCIF among other tumor shapes with the same volume. 
Equation (17) is able to legitimize this behavior of tumor 
shapes. The source term (the last term in the right hand 
side of Equation (17)) includes diffusion and convection 
terms. The convection term depends on the differences 
between IFP and P e ff based on starling s law. This pressure 
difference is close to zero (Figures 2 and 3) for all tumor 
shapes except for the prolate one and therefore the con- 
vection term only for prolate shape has non-zero value. 
Because of the non-zero values of Peclet number for pro- 
late shape in tumor region (Figures 11 and 12), the con- 
centration for this shape is affected not only by diffusion 
rate from vessels but also by convection rate from vessels. 
The convection rate leads to higher level of DCIF in prolate 
shape than that of other tumor shapes (Figures 6, 7, 8, 9 
and 10). The non-zero values of Peclet number for prolate 
shape are seen in an image based work of Zhao et al. as well 
[13]. In normal tissue, Peclet number shows that drug de- 
livery from microvessel to tissue is done by both mecha- 
nisms of transfer, convection and diffusion. 

The other effect of uniform pressure is zero pressure 
gradients in all tumor shapes except for prolate one. Be- 
cause of zero IFP gradient and based on Darcy s equation, 
IFV is close to zero in tumor tissue except for prolate 



shape. The close to zero value of IFV is predicted in a few 
numerical studies such as Welter and Rieger [36]; Roy and 
Riahi [37]; and experimental results of Hompland et al. 
[35]. Zero IFV results in a negligible convection transport 
(the second term of the right hand side of Equation (17)) 
and consequently the convection transport does not affect 
drug distribution; and the diffusion transport (the first 
term of the right hand side of Equation (17)) is the only 
reason of drug transport in all tumor shapes except for 
prolate one. Therefore, non-zero IFV in prolate shape, also 
seen in Zhao et al. [13], is another reason of higher DCIF 
values in this tumor shape with respect to other shapes. 

The sharp pressure gradient (Figures 2 and 3) and 
highest value of IFV (Figures 4 and 5) in the inner 
boundary for all tumor shapes increases drug transport 
and make a bump at this boundary for DCIF. 

The sensitivity analysis of effective pressure shows that 
the effective pressure does not have too much effect on 
DCIF. In all tumor shapes, DCIF for a wide range of ef- 
fective pressure changes smoothly; however, in prolate 
shape, this change is sharper and increases by effective 
pressure (Figure 13c and d). As mentioned, in the tumor 
region for all tumor shapes except for prolate one the con- 
vection rate, which depends on P e p has a negligible role 
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Volume (cm 3 ) 



Volume (cm 3 ) 



Figure 14 The effect of different values of tumor size on IFP, IFV and DCIF. a) The variation of maximum IFP in tumor region, b) The 
variation of averaged IFV at tumor boundary, c) Average of DCIF in tumor region for continuous injection, d) Average of DCIF in tumor region for 
bolus injection. 



on drug distribution; therefore, P e ff cannot have significant 
effect on drug concentration. However, in prolate shape, 
increasing P e ff increases pressure difference between IFP 
and P e jf%.no\ consequently the convection rate from vessels 
in the tumor region; therefore, DCIF in prolate shape is 
sensitive to P eff changes. 

The tumor volume shows more effects on IFP, IFV, and 
DCIF than other investigated parameters such as effective 
pressure and hydraulic conductivity. The increasing tumor 
volume increases significantly IFP (Figure 14a). The de- 
pendency of IFP to tumor volume is observed in experi- 
mental study of Gutmann et al. [38], Hompland et al. [39], 
and Leguerney et al. [40], as well. When the tumor vol- 
ume is in the order of 1 cm 3 , the sensitivity of IFP to 
tumor size decreases. The independency of IFP to tumor 
volume in in large tumors is observed in the study of 
Leguerney et al. [40], as well. In their work, IFP changes 
very slowly with tumor volume. Since the high IFP is in- 
troduced as the main barrier of drug delivery to tumors, 
IFP increasing with tumor volume leads to DCIF decrease 
in these tumors. This reduction of DCIF with tumor size 
is observed in Au et al. [41], as well. Since IFP reaches 
around the effective pressure with increasing tumor vol- 
ume, the convection rate is vanished and the diffusion rate 



reaches a constant value, and consequently the sensitivity 
of DCIF to tumor size reduces. Lower IFP in the small 
tumor sizes leads to increase the convection rate of source 
term in Equation (17). Therefore, it is expected to have a 
better drug distribution in small tumors. 

Results show that IFV, IFP, and DCIF are sensitive to 
tissue hydraulic conductivity changes. The hydraulic 
conductivity is appeared only in Darcys law. This par- 
ameter has a direct effect on IFF and indirect effect on 
DCIF. Theoretical analysis shows that the ^ in tumor 

region is proportional to 1-^/sinh^^ [8] (i< is hy- 
draulic conductivity, see material section). In low values 
of /c, the ^ / sinh is negligible and P t is close to P e ff 

(Figure 15a). Increasing /c, increases the ^/sinh^^ 

and leads to sharp decrease in IFP. This dependency is 
also observed by McCarty and Johnson [42], For high 
values of /c, IFP is very low and negligible in comparison 
to P e ff. IFP reduction from effective pressure increases 
IFV around 5 times (Figure 15b). The hydraulic conduct- 
ivity affects DCIF by convection rate from vessels (as 
mentioned this value depends on difference between IFP 
and P e ff)' In low values of /c, since IFP is equal to P e p the 



Sefidgar et al. Journal of Biological Engineering 2014, 8:12 
http://www.jbioleng.org/content/871 /1 2 



Page 9 of 13 




effect of convection rate is not significant and DCIF re- 
mains constant. Increasing k increases the convention 
rate and consequently DCIF. 

When hydraulic conductivity increases two to three or- 
ders of magnitude, the mean values of DCIF are two times 
greater than the average of DCIF for baseline values in 
Table 1. However, after a specific value of hydraulic con- 
ductivity, DCIF changes smoothly and reaches a constant 
value because IFP is very low and convection rate only de- 
pends on P e ff. 

Conclusions 

A numerical approach which couples the mathematical 
model of the lymphatic system and the interstitial flow 
with the mathematical model of solute transport demon- 
strates that DCIF is affected by two transport mecha- 
nisms, convection and diffusion. 

Drug convection and drug transport from microvessel 
to tumor are blocked by high interstitial pressure (IFP) 
which is uniformly distributed in most part of the 
tumor. The large pressure gradient results in an out- 
ward convective flow that washes out the drug extrava- 
sated from microvessels at the tumor periphery. This 



study shows that when there is IFP gradient in the 
tumor region instead of uniform IFP distribution which 
occurs in some tumor shapes, DCIF is greater than that 
of the uniform one. 

As the effects of osmotic and intravascular pressure 
on convection rate are negligible in most of tumor 
shapes, the dependency of DCIF to these parameters is 
very low. 

The hydraulic conductivity which is another considered 
parameter in sensitivity analysis has significant effect on 
drug distribution since it increases the convection rate 
from vessels. 

Method 

Mathematical model of interstitial fluid transport 

This section introduces the mathematical model of 
interstitial fluid transport in macroscopic scale [2,16]. 
Since normal and tumor tissue have characteristics the 
same as porous media, fluid flow behavior is defined by 
coupling the fluid flow governing equations. The mass 
balance or continuity equation for steady state incom- 
pressible fluid in the porous media with source and 
sink of mass is [16]: 
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V.V; = (p B -(p L 



(3) 



V,; = -KVPi 



(8) 



where 

V{. The interstitial fluid velocity, 

(p B : The source term, extravasation from micro vessels, 
and 

(f) L : The drainage term, elimination by lymphatic system. 
In biological tissues, the fluid source is evaluated through 
Starling's law as follows [16]: 



<Pb 



L P S 



{P B -Pi-o s {n b -jii)) 



(4) 



where 

Pi : Interstitial fluid pressure, 
P B : Blood pressure in microvessel, 
y : The surface area per unit volume of tissue for 
transport in the interstitium, 
n B : Microvessel oncotic pressure, 
iii : Interstitial oncotic pressure, 

L p : The hydraulic conductivity of the microvessel wall, 
and 

a s : Osmotic reflection coefficient. 

and the lymphatic system is related to the pressure dif- 
ference between the interstitium and the lymphatic ves- 
sels and is considered only for normal tissues [43] : 



Mr) 



LplSl 
V 

0 



(Pi -Pi ) Normal tissue 
Tumor tissue 



(5) 



where 

(f) L : The volumetric flow rate into the lymphatic, 

The lymphatic filtration coefficient, and 
P L : The hydrostatic pressure of the lymphatic. 
The momentum balance equation in its general form 
can be written as Equation (6) [44] : 



Pl^7+(v*-V)v, 



-Pi + ft Vv/ + (Vv/) 



2{t 
' 3 
h 

K, 



V.v< 

0*+ 



(6) 



where 

K: The permeability of the porous medium, 

p : The density, 

\i: The viscosity, and 

F: The volume forces. 

Since interstitial fluid is a Newtonian fluid and has low 
velocity through the tissues, Equation (6) in the steady 
state condition is simplified to Darcys law [16]: 



VP; 



(7) 



The K/u is defined as hydraulic conductivity of the in- 
terstitium, k: 



Combination of momentum equation (Equation (8)) 
and the continuity equation (Equation (3)) results in 



-v.(KVPi) = 4> B -4> L 



(9) 



When k is constant, the interstitial pressure can be 
expressed by 

-kV 2 Pi = 4> B -<p L (10) 
By substituting Equations (4) and (5) into Equation (10): 



(Pb -Pi-<?s fab -Hi) ) - ^ P y 1 [Pi-Pl) Normal tissue 



{P B -Pi-G s (n b -7ii)) 



Tumor tissue 

(ii) 



Macroscopic solute transport 

Molecular transport in tumors is based on the conserva- 
tion laws for mass and momentum. The interstitial trans- 
port of drug is governed by the convection diffusion 
equation; therefore, the general equation for the mass bal- 
ance of solutes can be written as: 



3C 



V.(D^.VC)-V.(v/C) + (Qb-Ql) 



(12) 



where 

C : The solute concentration based on tissue volume, 

®£ : The rate of solute transport per unit volume from 
microvessel into the interstitial space, 

<& L : The rate of solute transport per unit volume from 
the interstitial space into lymphatic vessels, and 

D e ff : The effective diffusion tensor. 

For an isotropic and uniform diffusion in tissues, equation 
(12) can be written as: 



dC 



■ D eff v 2 c-v-(v f c) + (a> B -o L ) 



(13) 



The solute transport rate across the lymphatic vessels 
can be considered as [15] 



*-{*o 



C Normal Tissue 
0 Tumor Tissue 



(14) 



The solute transport rate across the microvessel can 
be represented by Patlak equation [45]: 



PS Pe 
O b = <p B (l-a f )C P + — (C P -C)^- i 



(15) 



where 
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Pe 



PS 



(16) 



(p B : The fluid flow rate per unit volume of tissue 

across the microvessel wall, 
Of : The filtration reflection coefficient, 
P : The microvessel permeability coefficient, 
S/V : The microvessel surface area per unit volume of 

tissue, 

C p : Solute concentration in the plasma. 
By substituting Equations (15) and (14) in to 
Equation (13): 



Generating geometry 
and meshing 




Implementing 
governing equation of 
IFF on domain 










Calculating IFY and 
IFP in domain 




Using EB-FV 
numerical method 












Implementing 
governing equation of 
solute transport 




Calculating DCIF in 
domain 


1 


I 



Figure 17 Algorithm of numerical simulation used for 
calculating interstitial fluid flow parameter (IFV, IFP) and solute 
transport parameter (DCIF). 



dC 
dt 



= A # V 2 C-V-(v/C) 



PS 



C \ Pe 



C Normal region 



= D e/ V 2 C-V-(v/C) 



+ 1 <Pb (1-0/) Cp + y I c P - J I Tumor re s ion 



(17) 

Boundary conditions 

A tumor surrounded by normal tissue is considered in this 
study. The tumor shape is considered ellipsoid. A 2D sec- 
tion of the geometry and boundaries are shown in Figure 1. 
The three boundaries are indicated for geometry: a) the 
center of tumor, b) the boundary between tumor and nor- 
mal tissue, named inner boundary, c) the boundary at the 
outer edge of geometry, named outer boundary. The 
appropriate boundary conditions are implemented for 
Equations (11) and (17). The no flux boundary condition 
is applied at the center of the tumor; i.e. [16], 



VP, = 0forr = 0 



D € VC + ViC = 0forr = 0 



(18) 



(19) 



The continuity of pressure and velocity for Equation 
(11) and concentration and its flux for Equation (17) are 
considered as appropriate boundary conditions for inner 
boundary: 



-KtVPi 



Or — 
Pi\^ 



-K„VPi\ 



(20) 



Sefidgar et al. Journal of Biological Engineering 2014, 8:12 
http://www.jbioleng.org/content/871 /1 2 



Page 12 of 13 



{p{ ff VC + vtC) \a- = (p n eff VC + vtC) \ n - (21) 
c \n~ = c \n + 

where D~ and >Q + indicate the tumor and normal tissue 
at the inner boundary. 

For outer boundary, that the interstitial pressure is 
constant; the Dirichlet type of boundary condition must 
be applied [16]: 

Pi = for outer region (22) 

And for concentration, the open boundary condition is 
used in the outer region [46]. The Open Boundary is 
used to set up mass transport across boundaries where 
both convective inflow and outflow can occur and de- 
fined by Equation (23): 

-n-VC = 0 (23) 

where n is the normal vector. 

Geometry 

To investigate the effect of tumor shape on drug deliv- 
ery, 5 different shapes are considered. The 3D funda- 
mental shape of tumors is assumed to be an ellipsoid in 
different studies. The assumption of considering ellips- 
oid tumor shape is investigated in many research such 
as breast cancer [47], prostate cancer [48,49] cervical 
cancer [50], glioma cancer [51], and others [52]. Based 
on this mentioned reason, ellipsoidal shapes of tumors 
are considered in this study. Different shapes are pro- 
duced by changing ratio of two radii of ellipsoid shown 
in Figure 1. This ratio (A = bid) is changed from 0.1 
(prolate) to 10 (oblate). In all shapes shown in Figure 16 
the volume of the tumor is remained constant and equal 
to the volume of spherical tumor with radius R. The 
baseline value of R is 1 cm. The values of R used in sen- 
sitivity analysis are changed from 0.1 cm to 2.5 cm. This 
range is obtained from the literature and is close to the 
experimental observations [40]. 

Model parameterization 

The interstitial transport properties for normal and 
tumor tissue are listed in Table 1. These values are used 
as baseline and some of them are investigated and chan- 
ged in specified ranges for sensitivity analysis. 

The parameters of solute transport model taken from 
Baxter and Jain [8] are listed in Table 2. Although, the nu- 
merical model is applicable for any type of solute, in 
present study the properties of Fragment antigen-binding 
(F(ab')2) as a sample is used. 

Numerical solution 

A systematic flow chart is illustrated in Figure 17 to clar- 
ify the computational techniques involved in this paper. 



The criterion for the convergence of iterative solution 
based on element-based finite volume (EB-FV) method 
is to reduce the residual by 6 orders of magnitudes. The 
details of numerical method is mentioned in our previ- 
ous works [2,16,53]. 
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